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J.— V ', Abstract 

^^ , We study the settUng of sohd particles within a viscous incompressible fluid contained in a two-dimensional 

^ ' channel, where the mass density of the particles is slightly greater than that of the fluid. The fluid-structure 

^P( , interaction problem is simulated numerically using the immersed boundary method, with an added mass 

^^ ' term that is incorporated using a Boussinesq approximation. Simulations are performed with a single circular 

04 , particle, and also with two particles in various initial configurations. The terminal settling velocities for the 

particles correspond closely with both theoretical and experimental results, and the single-particle dynamics 

r^ ' , reproduce expected behavior qualitatively. The two-particle simulations exhibit drafting- kissing-tumbling 

^~^' dynamics that is similar to what is observed in other experimental and numerical studies. 

I ^ Keywords: immersed boundary method, particle suspension, sedimentation, settling velocity, 

3 ' fluid-structure interaction 

2010 MSC: 76T20, 65M08 
C/3 ' 
O 



1. Introduction. 

Particulate flows involve a dynamically evolving fluid that interacts with solid suspended particles, and 
arise in a wide range of applications in natural and industrial processes [9|] . We are particularly interested in 
►^ the gravitational settling or sedimentation problem, in which the suspended solid particles have large enough 

.,—1. , mass that they settle under their own weight. Sedimentation is observed in many applications, including 

(^ ■ flow of pollutants in rivers and the atmosphere, tea leaves settling to the bottom of a teacup, industrial 

00 , crystal precipitation, mineral ore processing, and hail formation in thunderclouds, to name just a few. 

^^ ' There is an extensive literature on experimental, theoretical and computational studies of particulate 

^^ , flows involving sedimentation. We make no attempt here to perform a comprehensive review, but will 

^D ' rather highlight a few of the more important results. Experimental studies of sedimentation have had a l ong 

^^ \ history including the earlier work of Richardson and Zaki |40| and extending to more recent years ISl. Il7. 

1271 . |28| . Many analytical and approximate solutions have been developed to explain the behavior of settling 



suspensions, especially in the dilute limit where there are only a small number of particles. Back in 1851, 
k>( ' Stokes [43| derived an analytical solution for a single particle settling within an unbounded fluid, and many 

^ \ other authors have since extended these results to other more practical sedimentation problems J6l.ll2l.l2ll.l47|. 

Cu ■ More recently, many numerical approaches have been applied to simulate settling particles, including the 

finite element method 16|, [l9|, |26|, |33[ , lattice-Boltzmann method [ISJ, |30|, |39| , and boundary element method 
24 136| . The underlying feature of these numerical methods is that the fluid flow is governed by the Navier- 
Stokes equations whereas the particles are governed by Newton's equations of motion. The hydrodynamic 
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forces between the particle and fluid are obtained from the solution of this coupled system, which typically 
requires either complex interfacial matching conditions at the fluid-particle interface, or else some form of 
dynamic boundary-fitted meshing. In any case, these methods tend to be complex and extremely CPU- 
intensive, especially for three-dimensional flows. 

One numerical approach that has proven to be especially effective for solving complex fluid-structure 
interaction problems involving dynamic moving structures is the immersed boundary (or IB) method. This 
approach has been used extensively to simulate deformable structures arising in problems in biofluid me- 
chanics 35] . Wang and Layton [48[ have recently used the IB method to simulate sedimentation of multiple 
rigid ID fibers suspended in a viscous incompressible fluid, and several other authors have applied the the 
IB approach to solve related sedimentation problems [7|, |l3|, |25|, |46|, |49j . 



The IB method is a mixed Eulerian-Lagrangian approach, in which the fluid equations are solved on 
an equally-spaced rectangular mesh, while the moving solid boundaries are approximated at a set of points 
that moves relative to the underlying fluid grid. In the original IB method, the effect of these immersed 
boundaries is represented as a singular force that is computed from the IB configuration and which is then 
spread onto fluid grid points by means of a regularized delta function. The added mass due to a sedimenting 
particle can also be distributed onto the fluid in a similar manner. With the exception of the papers by 
Wang and Layton 48* and Hopkins and Fauci [25], the other authors mentioned above have employed a 
modification of this IB approach known as the "direct forcing IB method," wherein the force is an artificial 
quantity that is calculated directly from the governing equations so as to satisfy the velocity boundary 



conditions exactly on the immersed boundary (see [32] for more details) 



Our aim in this paper is to apply the original IB method to solving sedimentation problems, rather 
than the direct forcing approach. We restrict ourselves to a two-dimensional geometry, in which one or two 
particles with a circular cross-section settle under the influence of gravity within a rectangular channel that 
has vertical bounding walls. Although the IB approach has been applied to solve certain sedimentation 
problems, there has not yet been an extensive comparison to other results in the literature. Our pri mary 
aim is therefore to perform such a comparison to a number of experimental 4J, l50| , theoretical J14l . |45| , 



and numerical [l(x\ studies, in order to ascertain the validity of the IB approach in simulating sedimentation 
problems. Although we focus here on solid particles, the long-term goal of our work is to develop a numerical 
framework that can be used to investigate the settling of highly deformable particles. 

We begin in Section[2]by describing the IB method and defining the forces used to simulate the presence of 
both settling particles and channel walls. Section[3]contains a review of previous analytical and experimental 
results on the settling velocity for a single particle in both unbounded and wall-bounded domains. We then 
perform a series of numerical simulations of sedimentation at small to moderate Reynolds numbers, and 
report the results in Sections S] and [Sj Most of the results appearing in this article are contained in the PhD 
thesis of the first author [18| . 

2. Immersed boundary method. 

The immersed boundary method is both a mathematical formulation and a numerical scheme. We 
begin in this section by describing the model equations that underlie the IB formulation for fluid-structure 
interaction. Following that, we discretize the equations and describe the numerical algorithm used to 
determine an approximate solution. Finally, we provide details on the specification of the discrete IB force 
density representing the channel walls and sedimenting particles. 

2.1. Model formulation. 

In this section we describe a two-dimensional IB model that is capable of capturing solid (and potentially 
deformable) elastic bodies with general shape and that move within a surrounding incompressible, Newtonian 
fluid under the action of gravitational force. The details of the IB force density used to handle a solid circular 
object in the presence of two parallel bounding walls are left for section [^31 All variables and parameters 
in this paper are stated in CGS units, unless otherwise indicated. 

Suppose that a moving elastic solid body F is contained within a fluid domain J7 as pictured in Figure [TJ 
In general, F may consist of several disconnected components, F — IJjFi, where each F^ can be a one- 




a 



Figure 1: A general immersed boundary configuration F : 
within a doubly-periodic fluid domain f2. 



Ui=i ^i consisting of several disconnected components immersed 



dimensional elastic membrane (parameterized by a single real parameter s) or an elastic solid region (whose 
specification requires two parameters, r and s). We denote the location or "configuration" of the immersed 
boundary by X(q,t) [cm], where q is a dimensionless IB parameterization that is used to represent either 
a scalar s or a vector (r, s), depending on the context. For simplicity, we assume that fl = [0,Lx] x [0,iy] 
is rectangular in shape and that periodic boundary conditions are applied in both the x- and y-directions. 
The effect of the elastic body on the fluid is to impose a force fj^ [g/cm^ s^] onto the adjacent fiuid 
particle at location x = X{q,t), which is incorporated into the incompressible Navier-Stokes equations as 
follows: 



V • M = 0. 



(1) 

(2) 



Here, u{x,t) is the fluid velocity [cm/s], p{x,t) is the pressure [g/cms^], x — {x,y) are the Eulerian 
coordinates [cm] for the fluid domain il, p is density [g/cm'^] and /i is dynamic viscosity [g/cms]. The IB 
forcing term in the momentum equations ((ij is represented by a force density FiB{q,t) [g/s^] that is spread 
onto the surrounding fiuid by means of a delta-function convolution 



fjB{x,t)= J FiBiq,t)S{x-X{q,t))dq, 



(3) 



where S{x) — 5{x)S{y) is the Cartesian product of two one-dimensional Dirac delta functions. 

Most papers in the immersed boundary literature assume that F has the same constant density pf as 
the surrounding fluid, and hence F is neutrally buoyant. However, for the particle sedimentation application 
considered here, we must take F (or at least portions of it) to have density ps > pj that is greater than that 
of the fluid. Consequently, the density of the fluid-solid composite material p[x,t) is a variable quantity 
that may also be written in terms of a delta function convolution as [52| 



p{x,t) = Pf +Ap{x,t) 



where 



Ap{x, i)= I M{q) 5{x - X{q, t)) dq. 



The quantity M{q) ^ is the added Lagrangian mass density due to F, with M = only for those 
components that are neutraUy buoyant. 

In all examples in this paper, we will take M = Mq (constant), and we also assume that the solid 
density is close to that of the fluid so that Ap <C p/ . Consequently, it is reasonable to apply a Boussinesq 
approximation as in [25| so that the extra intertial term involving Ap is neglected and the density on the 
left hand side of the momentum equations ([1]) is taken equal to the constant pf : 

CJU 

Pf-Q^ + Pf^ ■ ^^ = P"^^^ - Vp + fjB + fa- (4) 

The extra forcing term /^ derives from the force of gravity acting on the immersed boundary and can be 
written as [25[ 

fa{x,t)^-gkAp^-gkjM{q)S{x-X{q,t))dq, (5) 

where g — 980 cm/s^ is the gravitational acceleration and k = (0,1) is the unit vector in the vertical 
direction. 

Finally, the immersed boundary is assumed to move with the fluid so that 

f)X f 

— = y u{x, t) Six ~ X(q, t)) dx, (6) 

which is simply the "no-slip" condition for fluid particles located adjacent to the immersed boundary. 

In summary, the governing equations consist of (0), (JH)-®, with the IB force density being the only 
component that remains to be specified. Since it is easiest to write fjg in discrete form, we will first derive 
the discretized governing equations, after which we will provide a specification for the IB force. 

2.2. Numerical algorithm. 

The algorithm we describe next is a semi-implicit scheme that is closely related to the method outlined 
in [4l|. The fluid domain il is divided into an equally-spaced grid of points denoted by Xij = {xi,yj) — 
[ihx, jhy), with hx = L^/N^, hy = Ly/Ny, i = 1, 2, . . . , N^, and j = 1, 2, . . . , Ny. We consider a time interval 
[0, T] divided into equally-spaced points i„ — nAt with time step At — T/Nt and n — 0, 1, 2, . . . , A't. We 
may then define discrete approximations of the velocity and pressure m" and p" at points {xi, yj, tn). The 
immersed boundary F is similarly discretized at points Xi for (. — 1,2, . . . , Nf,, and the IB configuration and 
force density are approximated by X^ and F" respectively. 

Using the above notation, we introduce finite difference operators that approximate the spatial derivatives 
appearing in the governing equations. In particular, we define two one-sided difference approximations of 
the ^-derivative of a grid quantity Wij 

i^>,, = ^^^-ip^ and i^-^., ^ "-- 7'-^^ (7) 



as well as the centered approximation 



D>,, ^ "'^^'^^7-^'^ (8) 



Analogous definitions apply for the y-derivative approximations D^ , D^ and D^, and the gradient is 
replaced by the centered approximation V/j = (D^, -Dy)- Finally, the delta function appearing in the integral 
terms is replaced by the regularized function 



i-xii-y 



^1 + cosf^n, if |r| ^ 2, 



where 

{r)^{AV '-\2JJ' -'■'-;' (10) 

0, otherwise. 

We are now prepared to state the immersed boundary algorithm. In any given time step, we assume that 
values of the velocity u^~^ and IB configuration X"~^ are known from the previous step. These quantities 
are evolved to time t„ using the following procedure: 

1. Compute the IB force density F"g^ based on the configuration X"~^ as described in section [^751 

2. Spread the IB force to the fluid grid points using a discretization of the integral in (|3]) 

e=i 

and a similar approximation of the integral in ([5]) yields a formula for /^^ . The scaling factor Ai, in 
both cases is inversely proportional to the number of IB points (Nf,) and has a different interpretation 
depending on whether the immersed boundary is a ID fiber (channel wall) or a 2D solid block (circular 
particle). In the case of a fiber Ai, is a length, while for a solid region Ai, is an area; in both cases, the 
factor Ab ensures that the formula dill) scales properly with the number of IB points and that it is a 
consistent approximation of the corresponding integral. More details on the precise form of (fTTj) and 
the specification of At are provided in section 12.31 

3. Integrate the incompressible Navier-Stokes equations using a split-step projection scheme: 

(a) Compute an intermediate velocity u] ■ by applying the elastic and gravitational forces on the 
immersed boundary: 

/..(I) ..n-l\ 

IB.i.j T" JG.i,] (12) 

(b) Apply an ADI discretization of the advection and diffusion terms: 

Pf (%^ + <j'D^<l) - PDtD-^, (13) 

/ (3) __ (2) \ 

Pf [^'^'^+<j'^<i) -PD-D^^. (14) 

These equations represent a sequence of tridiagonal solves for u\ and ui . 

(c) Project the intermediate velocity u\ onto the space of divergence-free vector fields by: 
i. Solving the pressure Poisson equation 

V„-V,.p.,, = ^V,.«g. (15) 

Note that V/,. • V/i represents a wide finite difference stencil for the Laplacian involving 
the pressure values pij, Pi-2,j, Pi+2,j, Pi,j-2 and Pi,j+2- Owing to the periodic boundary 
conditions on il, this equation is solved most easily and efficiently by means of the discrete 
Fourier transform, and making use of the FFT algorithm. 




ii. Updating the velocity according to 



4. Evolve the immersed boundary using 



,(3) 



JJ 



^Y^, 



At 
Pf 



■ ^hPi. 



M2_^ulj5h{xij 



X g ) hxhy. 



(16) 



(17) 



This simple semi-implicit time discretization described above yields a solution that is first-order accurate 
in time. And although all spatial derivatives are approximated using second-order finite differences, the 
method is still only first-order accurate in space owing to errors in velocity interpolation near the immersed 
boundary that arise from the use of the regularized delta function. It is straightforward to increase the 
temporal accuracy to second order using an algorithm such as that proposed by Lai and Peskin 31], but it 
is much more difficult to increase the spatial accuracy [20[ . Since the focus of the current study is to validate 
the general IB approach in the study of particle sedimentation, we have chosen to employ the simple scheme 
above, and leave for future work the implementation of higher order extensions to the algorithm. 

2.3. Discrete IB force density for particle and channel walls. 

We begin by describing the geometry for the particle sedimentation problem. Referring to Figure [21 we 
take a rectangular fluid domain of size L^ x Ly and place two vertical immersed fibers representing the 
channel walls a distance W < L^ apart, symmetric relative to the channel centerline, and separated from 
the domain boundary by a narrow strip of fiuid. With periodic boundary conditions applied on all sides of 
the domain, the channel walls naturally connect to each other across the top and bottom boundaries. A 
single, solid, circular particle of diameter D is initially located at the center of the channel. Later on, we 
will consider other initial configurations with one and two particles, but for now this will suffice to illustrate 
the calculation of the IB force density. This circular particle in 2D may be thought of as corresponding in 
3D to a cross-section of a solid cylinder with infinite length. 




Figure 2: Initial geometry for the gravitational settling problem. The parallel channel walls are denoted by dashed vertical 
lines, separated by a distance W. The first test case has a single solid particle of diameter D that is released at time t = 
along the channel centerline (which is indicated by a dotted line). 

In our sedimentation model, the IB force density Fjb is the sum of two terms, Fjb — F^ + F'^, where 
F^ represents the force density generated by the channel walls and F"^ is that generated by the circular 
particle. These forces are discussed separately in the next two sections. 



2.3.1. Elastic force from the channel walls, F"^ . 

The vertical walls are discretized using an equally-spaced array of IB points that are initially at locations 
X;'"^ = (i(L^ - W), £h^) for the left wall, and X^'^ = (i(L^ + W), £h^) for the right wall, where the 
wall point spacing is /i^, = Ly/N^ and i — 1,2,..., N^. Each IB point is connected to a fixed "tether point" 
(at the same initial location) by a very stiff spring that exerts a force of the form 

F7'^-a^(X™'^-X,), (18) 

where ct^, [g/cms^] is the spring stiffness and X^ is the moving IB point location. Any motion of the IB 
point away from corresponding the tether point location generates a spring force that drives it back towards 
the target, so that as long as cr^ is chosen large enough the wall points can be made to mimic a rigid 
structure. We emphasize that tether points neither move with the fluid nor generate any force themselves. 
A similar expression is developed for the force density at the right wall points, i^^' so that the total wall 
force density may be written as 

F™=^(Ff^ + F,"''^). (19) 

1=1 

The natural choice of scaling factor in the force spreading step (jlip is the wall point spacing, Af, — hw 

2.3.2. Elastic force from the particle, F^ . 

The circular particle is represented by a collection of N^ Lagrangian points that lie on its circumference 
and throughout its interior. We make use of the unstructured triangular mesh generator DistMesh [3^ that 
generates a nearly uniform triangulation such as that shown in Figure |3l The nodes of the triangulation 
are the IB points Xi, for € = 1, 2, . . . , A^c, while the edges of the triangles define a network of springs that 
maintains the shape of the particle. In addition to bearing IB spring forces, the network nodes are also 
employed in equation ([5]) to distribute added mass throughout the particle. In practice, we generate the 




Figure 3: Uniform triangular mesh generated by dlstmesh2d. 

triangulation by calling the Matlab function distmesh2d with the "scaled edge length function" hunif orm 
(a function provided by the authors that attempts to find a mesh that is as uniform as possible). We also 
set the "initial edge length" parameter equal to ^ xmx\{hx, hy), which ensures that the mesh obeys 

max \Xk - Xe\ < - iam{hx, hy), 



which is a standard "rule of thumb" that avoids leakage of fluid between IB points 



This form of particle discretization should be compared with the more common IB approach that uses 



an open circular ring of points with a freely-moving fluid inside, such as in [311 . |49|. This approach has 
been criticised [22| for generating non-physical fluid motions inside the particle and in some cases leading 
to significant deviations in the shape of the particle. In contrast, our discretization of the particle interior 
with a network of IB springs suppresses this spurious fluid motion and also helps to maintain the rigidity of 
the particle boundary. 

We now define the spring forces that act on the network, following the development of Alpkvist and 
Klapper for viscoelastic biofilm structures [2|. Let di^rnif) — ^lit) — ^mit) be the vector joining two IB 
points labeled H and m, and let df,m(i) — |d^,m(i)| be the corresponding distance. We assume that the 
spring network is initially in equilibrium (i.e., zero force) so that all springs have a resting length equal to 
their initial length, di^m{Q)- Let I be an incidence matrix whose entries li^m are either I or depending 
on whether or not points i and m are connected, respectively. Then the force density acting on the £*'* IB 
point in the network is 

F'i^a, J2 hm-r^{df.,,niQ)~de^m), (20) 

=1 

#0 



m— 1 



where the sum is taken only over those m for which Xm is connected to Xi in the network. We have also 
assumed that the spring stiffness CTc [g/cms^] is constant for all network connections. The total elastic force 
density generated by all IB points making up the circular particle is then given by 

F^ = Y,Fl (21) 

1=1 

The appropriate scaling factor for the force integral ([TT|) is the average area of a triangular mesh cell, 
Ab = 'k{D/2)'^/Nc. a similar approach was employed by Hopkins and Fauci [25J to simulate a suspension of 
microbial cells that they treated as point particles. 

3. Approximate formulas for settling velocity. 

We next review some of the existing analytical and experimental results on the settling of a single particle 
falling under the action of gravity. The study of a spherical particle in an unbounded fluid medium in 3D 
is a classical problem that was considered by Stokes [43|, who obtained a formula for the settling velocity 
that is now known as Stokes' law. We will first state Stokes' result and then modify it for a circular particle 
in 2D, which corresponds to an idealized "infinite cylinder" in 3D. We then consider the case of a circular 
particle falling in a bounded fluid domain between two vertical walls and then review several of the most 
commonly-used formulas for the "wall-correction factors" that have been obtained from either fitting to 
experimental data or using approximate analytical techniques. A fairly extensive overview of settling for 
cylindrical particles, including many of the wall correction formulas reported in the literature, is given by 
Champmartin and Ambari [10|. 

3.1. Stokes' law for a spherical particle in 3D. 

There are two main forces acting upon a massive particle settling in a fluid: the gravitational force Fg, 
and the drag force Fd due to the "friction" between the particle and the fluid. A particle that is initially 
at rest will accelerate under the action of gravity, and as the particle begins to move through the fluid it 
experiences a drag force in the direction opposite to its motion that increases with the speed of the particle 
relative to the fluid. If the drag force becomes large enough that it equals the gravitational force, then the 
two forces are in balance and no further acceleration occurs. The particle velocity in this equilibrium state 
is known as the settling or terminal velocity. 



We take a sphere of diameter D and density ps whose added mass relative to the fluid is ^tt (-j) {pg — pf). 
The net gravitational force acting on the sphere is 

and the corresponding drag force is 



^s-yly) 9{Ps-Pf), (22) 



^x2 



F'i^2^dPf'^^y-^) ' (23) 

where Cd is the drag coefficient for a sphere and V is the velocity of the sphere relative to the fluid. The 
settling velocity Vg corresponds to the long-term steady state in which drag and gravity forces are in balance, 
so that Fd = Fg. By equating (1^^ and (P5|) . we can solve for 



keeping in mind that the drag coefficient on the right hand side also typically depends on the settling velocity, 
Vs. Indeed, we know from 3] that the drag coefficient for a sphere can be approximated for small Reynolds 
number by 



where we have taken 



based on the particle diameter. Substituting this expression into ([24]) and solving for Vg we obtain Stokes' 
law 



Re^Pl^^, (26) 



_ gD\ps-pf) 



which is valid for Re < 0.1. 



3.2. Settling velocity for a circular particle in 2D. 

A similar argument may be used to derive the corresponding expression for a circular particle in 2D. 
We begin by considering a cylinder with diameter D and length £ and take the limit as i? ^^ oo in order to 
obtain a result that is relevant to our 2D geometry. The immersed cylinder has an added mass (relative to 
the fluid) of m = TT (y) £{ps — pj) ior which the net gravitational force is 

Fg^^giD'{ps-pf), (28) 

and the drag force is 

Fd = ^CdPfV^Dl. (29) 

Notice that the cross-sectional area factor tt {D/2) for the sphere from (j23l) is replaced by D£ for the cylinder, 
and the tildes are used here to denote cylindrical quantities. We also make use of the drag coefficient for a 
cylinder from j3| 



Btt 
Reh[ 



Crf = „„,„ ,7.4V (30) 



Re 



Ps 


Ap 


v; 


1.01 


0.01 


0.0024 


1.02 


0.02 


0.0045 


1.03 


0.03 


0.0066 


1.04 


0.04 


0.0086 


1.05 


0.05 


0.0105 



Table 1: Settling velocities Vs for a cylindrical particle, obtained by solving equation II32II with pf = 1 and D = 0.018. 



which holds when i?e <C 1. 

The setthng velocity for the cylinder is then obtained by equating the gravitational and drag forces in 
dMl) and dSll), which yields 



V,^ 



hgPjps - Pf) 
2CdPf 



(31) 



Observe that the factor of length i cancels in the above expression, so that this same expression is valid also 
for the 2D geometry in the -^ — t- oo limit. Furthermore, this expression is the same as that for the sphere in 
f^M except that the factor •y/4/3 is replaced here with •^7r/2, and of course the cylinder drag coefficient is 
also different. When equations (I5ni) - (PT|) are taken together, they reduce to a nonlinear equation in Re 



f{Re) = Re 



PjgD^JPs - Pf) 



■"'S 



0, 



(32) 



which can alternatively be written as an equation in Vs- It is easy to show that the function f{x) is 
continuous on the interval < x < oo and has the following properties: 



lim f{x) 

r-S.O+ 



-oo, 



lim J{x) 



and f\x) > 0. 



Therefore, / is guaranteed to have a unique positive real root by the intermediate value theorem. 

Newton's method may be used to solve (|5^ for Re, and we find that any initial guess for Re suffices 
since the convergence is to rapid. Table [1] lists values of Vs from ([5^ for parameters D — 0.018 and ps 
ranging from 1.01 to 1.05. As expected, the settling velocity increases with particle density as in the Stokes 
case. 

3.3. Wall- corrected settling velocities. 

In this section, we summarize a number of formulas that approximate settling velocity for a particle in a 
bounded fluid domain that consists of a channel with two parallel, vertical walls separated by a distance W. 
It is well-known that the bounding walls exert an additional retarding effect on a sedimenting particle [J, [ll , 
14l . [23, 37, 45 , 50] so that the settling velocity is lower in a channel than in an unbounded domain under the 
same conditions. The effect of these wall interactions may be approximated by means of a "wall correction 
factor" A, that is usually expressed in the form [J, |5| 



X{k,Re) 



Fdjk) 



(33) 



where Fd{k) represents a drag force per unit length, Vc is the terminal settling velocity in the channel, and 
k — D/W is the dimensionless particle size with < fc < 1. 
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Note that the factor A(fc, Re) depends on both Reynolds number and particle size; however, it is well 
known !ll| that at either very low or very high values of 5R, A is nearly independent of the Reynolds number. 
In this study, we are concerned with the low Re regime and so it is reasonable to assume that A depends 
on k only. Various formulas for the wall correction factor have been reported in the literature (see I4|, for 
example) all of which have the same channel geometry as pictured in Figure [51 Some of the more common 
wall correction factors are listed below: 

• White [50| : carried out experiments with various wires and ebonite rods in a channel containing viscous 
liquids such as glycerin and paraffin. He obtained the following experimental fit for the drag force on 
a cylinder 

whose domain of validity is restricted to < fc < 0.2. 



• Faxen [IJ, |23|: derived an approximate analytical solution of the Stokes equations, from which he 
obtained 

— 47r 
^^'^^ " 0.9157 + In(fc) - 1.724A:2 + l.TSO/c^ - 2.406fc6 + 4.591fc8 ' ^^^' 

Some authors claim that this approximation is valid for k as large as 0.5 |4l |53|. However, others cite 



an upper bound of fc — 0.3 or even lower 37[ which is more in line with our numerical simulations (see 
Figure [5] in Section E]). 



Takaisi [45| : used an analytical solution of Oseen's equations to obtain the approximation 



'^'^ = 0.9156 tln(fc) ' (^«) 

which is restricted to < A: < 0.2. He also performed a comparison with White's experimental fit and 
showed that the two expressions match reasonably well when k < 0.05. 

If we now consider A(fc) to be a known function of the dimensionless particle size fc, then equation p3p 
can be solved for the drag force per unit length as 

Fd{k) = VcfiX{k). 

Equating this expression with the gravitational force 

Fg^lgD'iPs-Pf), 

we find the following formula for the confined (or wall-corrected) terminal settling velocity of a cylinder 

^^" 4M(fc) ■ ^ ^ 

In the next section, this expression will be compared with numerically simulated values for the three choices 
of A(fc) listed above. 

4. Numerical results: Single particle case. 

In this section, we concentrate on a single particle that settles under the influence of gravity. Two initial 
configurations are investigated: one a symmetric case in which the particle is released along the centerline, 
and a second asymmetric case where the particle is released from an off-center location. 
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Figure 4: A convergence study in the settling velocity Vb on a sequence of grids. The difference in values of Vs on two 
successively finer grids is plotted versus the grid spacing h. A straight line with slope 1 corresponding to a first-order method 
is shown for comparison purposes. 

We restrict ourselves to a low Reynolds number regime corresponding to Re < 7, where Re refers to 
a "final" Reynolds number that is based on the vertical velocity after a particle has achieved its terminal 
settling velocity. Unless otherwise indicated, we choose physical parameters pf = 1, ps — 1.01 and D — 0.08. 
The wall and particle IB spring stiffness values are taken large enough that the walls and particle boundary 
do not deform "too much" from their initial shapes - taking cr^ = Cc = 3 x 10^ keeps the relative error in 
these boundary shapes to within approximately 0.2%. 

Except for the convergence study in the next section, most of our simulations are performed at the same 
grid resolution of hx = hy = 0.0083 and with a time step of At — 10^^. We also select the number of IB 
points for the two channel walls (iY^,) such that the ratio of spacing between wall points to fluid grid size is 
hf^/ uim(hx, hy) ~ ^ - this ratio is well within the factor of ^ that is recommended to avoid leakage of fluid 



between IB points [35|. We also adjust the number of IB points for the particle (Nc) until the initial mesh 



computed by DistMesh satisfies the same criterion. 

4-1- Convergence study 

We begin by performing a convergence study that validates the spatial accuracy of our numerical method. 
As mentioned earlier in section 12.21 the IB algorithm being employed here is well-known to be first order 
accurate in space. To verify this result, we select a sequence of fluid grids with N^ = Ny = 56, 112, 224 and 
448 on a square domain with side length L^ = Ly = 1, and use the settling velocity Vs as a representative 
measure of the solution for each case. The difference between values of Vs on successive grids is calculated 
and the results are plotted in Figure HI which demonstrates that our numerical solution converges as the grid 
spacing is reduced. The curve is nearly a straight line on a log-log scale, and the slope of 0.79 obtained from 
a least squares flt suggests that our implementation of the IB method is close to the expected flrst-order 
accuracy. Similar convergence rates are observed for other quantities such as fluid velocity, IB position, etc. 



4-. 2. Comparison with Stokes' law. 

We aim next to validate the numerical method against the settling velocity Vs for a cylinder in an 
unbounded medium. However, we recall that our doubly periodic geometry implies that a single particle 
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L 


computed V 


1 


0.001230 


2 


0.001462 


3 


0.001565 


4 


0.001635 


5 


0.001673 


6 


0.001674 


7 


0.001674 



Table 2: Computed settling velocity as a function of domain size L for D = 0.018. 

actually corresponds to an infinite array of sedimenting particles. Therefore, in the absence of solid bound- 
aries or any other mechanism for dissipating energy, the net effect of gravity acting on such an infinite array 
of mass-bearing particles will be to accelerate the particles and the surrounding fluid indefinitely. This 
situation is clearly non-physical, and so instead we introduce walls into the domain that are situated "far 
enough" from the particle so as to minimize wall-particle interactions and yet still permit the particle to 
reach its natural terminal velocity. To this end, we take a square domain with side length L^ — Ly — L that 
contains two vertical walls separated by a distance W = L — 0.04, and perform a sequence of computations 
with successively larger L. In particular, we fix the particle diameter at D = 0.018 cm and vary L between 
1 and 7 cm. The fluid viscosity here and in the next section is fi = 1 g/cm s, which we remark is much larger 
than in later sections because a meaningful comparison to Vs is only possible at low Reynolds number. 

The computed values of settling velocity are summarized in Table [H from which we observe that as L 
increases V approaches a limiting value of roughly 0.001674 cm/s. The results have clearly converged on the 
largest domain size, but the limiting value of V is significantly less than the Stokes settling velocity from 
equation (PTj) . Vg = 0.00239 cm/s. We suspect that the discrepancy is due to a combination of effects arising 
from grid resolution, including the first-order dispersive errors in our numerical scheme and the increase in 
the effective thickness of the walls and particle owing to the delta function smoothing width (this last effect 
is discussed further in section l4.3.2[) . 

4-3. Single particle initially along the centerline. 

We next consider the channel domain pictured in Figure [2] wherein the particle is initially released along 
the center of the channel. In this case, the symmetry suggests that any forces generated by particle-wall 
interactions are balanced and so the particle should fall along the centerline without veering to either side. 
We perform a number of sensitivity studies that investigate the effect of parameters such as the fluid domain 
size [0,La:] X [0, Ly], density difference Ap, and relative particle size k on the settling velocity. As in the 
previous section, simulations are restricted to low Reynolds number by taking /.t = 1. 



4-3.1. Dependence of settling velocity on density difference Ap. 

For a fixed channel size with W = 0.98 and L^ = 1, we vary Ap between 0.01 and 0.17. The plot of 
settling velocity in Figure [S] shows that the wall-corrected settling velocity Vc increases linearly with Ap, 
which is consistent with equation (1571) . We also consider the effect of changes in the channel length by taking 
Ly G {3,10,16}, which indicates that the infiuence of periodic copies in the j/-direction is most greatest 
for the shortest channel {Ly — 3) which also exhibits the highest settling velocity. As the channel length 
is increased, the settling velocity decreases until by Ly = 16 the results appear to have converged and are 
incidentally closest to the result predicted by Faxen's formula ([55]) . 



4.3.2. Dependence of settling velocity on particle size k. 

As mentioned earlier in section 13.31 the work of Faxen, White, Takaisi, and others suggests that A 
depends only on k at low values of Reynolds number. This motivates our next sensitivity study of the effect 
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Figure 5: Settling velocity in a channel as a function of Ap. parameter values: Lx = 1, Ly £ {3,10,16}, D = 0.1, W = 
0.98, k = 0.102. 

of dimensionless particle size, for which we again fix the channel width at Ty = L^; — 0.04, and then vary k 
by choosing values of particle diameter D £ [0.018,0.96]. Alongside our computational results in Figure IHl 
we have displayed corresponding estimates of the confined settling velocity Vc calculated using equation 
([57| with the three wall correction factors ([M1) - ([5S|) . The unbounded cylindrical settling velocity Vs is also 
included for comparison purposes, which clearly diverges from the wall-corrected values away from fc = 0. 
For these simulations, Re was in the range [1.2 x 10"^, 1.8 x 10^^]. 

Our computed results match most closely with Faxen's formula, which is most often cited as the most 
accurate approximation for the wall-corrected settling velocity. We also performed a study of the effect of 
changes in the channel length Ly, in order to determine the effect of any possible interference between vertical 
periodic copies of the particle due to the periodicity assumption. Results for Ly G [1, 16] are displayed in 
Figure [7] which clearly show that our computed settling velocity converges to a value quite close to that 
predicted by Faxen's formula when k < 0.2. 

On the other hand, there remain significant deviations between Faxen's results and our computations 
for values of k larger than 0.2, as demonstrated in Figure [5] Faxen's settling velocity levels out and attains 
a maximum near k = 0.3 and then falls to zero near k = 0.6. In contrast, our computed settling velocity 
reaches a maximum that is roughly 60% larger {Vc ~ 0.035) and plateaus for k roughly in the range [0.4, 0.8]. 

Our computed settling velocity only drops to zero when k is very close to 1, which is easily justified 
since the particle must come to a stop as it come into direct contact with the stationary walls. However, 
our results in Figure [S] show that Vc actually tends to zero not at fc = 1 but rather k « 0.96. The reason for 
this apparent reduction of 0.04 in the channel width is that the approximate delta function in our numerical 
scheme has a finite smoothing width that has the effect of introducing an extra "effective thickness" to both 
the walls and the particle. The numerical simulations in [42l] show that when using the cosine delta function, 
the effective thickness of an immersed boundary is approximately l.Qh, where h is the fluid grid spacing]]. 
Consequently, a particle with diameter D should have an effective diameter of roughly D^ff Ki D + 3.2h, 
while the walls should each extend an additional distance of 3.2/i into the channel. Taken together this 



^Note that the effective thickness depends on the choice of regularized delta function. Bringley [g| computed an effective 
thickness closer to 1.25h for a different but closely-related approximate delta function. 
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Figure 6: Variation of settling velocity with k for the computed results and various analytical expressions. Parameter values: 
Lx = 1, Ly = 16, fi = 1. 
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Figure 7: Effect of increasing the channel length Ly and hence reducing the impact of the periodic copies in y. Parameter 
values; ^ = 1, Lj; = 1, D = 0.2, W = 0.98, k = 0.2041. 



15 



Faxen eighth order 
Our simulations 




Figure 8: Plot of settling velocity versus k, comparing our simulations to Faxen's approximation. Parameter values: L^ = 1, 
Ly = 16, A* = 0.01. 



suggests a total reduction of 6.4/i in the effective channel width, which for h 
0.053. This is not far away from the observed reduction of 0.04. 

We summarize the behavior from our numerical simulations as follows: 



0.0083 equals approximately 



• For small particle diameters corresponding to k G [0, 0.2], the particle is far enough from the channel 
walls that the retarding effects of wall drag are not as prominent. In this range, the dependence of the 
settling velocity is roughly proportional to fc, which is consistent with Faxen's result. 

• For intermediate values of fc, roughly in the range [0.4,0.8], the settling velocity has attained a max- 
imum value and remains approximately constant. For these particle sizes, the interactions with the 
walls are at long range and are mediated by the fluid. 

• For values of fc G [0.8, 1.0], the particle is very close to the walls, giving rise to close-range interactions 
that slow the particle significantly. 

Of course, the validity of Faxen's approximation is limited to k < 0.2 and so it is no surprise that our results 
differ so much for larger k. 

4-4- Single particle initially off-center. 

In this section, we consider an asymmetry initial condition in which the particle is released from an 
off-center location. In Figure [HI the initial configuration labeled "t = s" shows the particle a distance 
W/2 to the left of center. The diameter of the particle is taken to be D = 0.08 cm, the channel length is 
Ly = 3, and we consider two different channel widths, W = AD and 8D. We also vary Reynolds number 
by taking values of the viscosity ^ G [0.006,0.018] g/cms. This choice of parameters allows us to draw 
a comparison with the analytical and experimental results reported by Sucker and Brauer [44|, as well as 
numerical simulations of Feng et al. [16[. 

The settling dynamics are pictured in Figure [3] for W = AD and Re ~ 4.9. As the particle falls, it 
initially drifts to the right toward the channel centerline, eventually attaining its terminal settling velocity 
there. The plot of particle trajectories in Figure [10] shows that the particle actually undergoes a damped 
oscillation about the channel centerline with an initial overshoot. Simulations were also performed for three 
other Reynolds numbers. Re — 2.2, 3.7 and 4.4, and the corresponding trajectories in Figure [TUl show that 
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Figure 9; (Left) Settling dynamics for a single particle released off-center, in a channel of width W = AD and Re = 4.9. The 
solid line drawn through the center of the particle highlights the rotational motion. (Right) Velocity vector plot at time t ss 2 s. 



increasing Re leads to larger oscillations about the centerline. For the smallest value of Re = 2.2, the 
particle trajectory undergoes a nearly nionotonic approach toward the centerline; in an analogy with simple 
harmonic oscillation, this behavior can be described as an overdamped oscillation. 

In addition to the vertical and horizontal translations of the center of mass, the particle also undergoes 
a small-amplitude rotational motion as it settles, which can be seen by tracking the progress of the straight 
line drawn through the center of the particle in Figure [SI This rotation can be more easily seen in the 
plot of angular velocity in Figure [TT] for the Re = 4.9 case. Initially, as the particle drifts from its starting 
location toward the centerline, it experiences a slight counter-clockwise rotation. As the particle approaches 
its equilibrium horizontal location, the rotation slows and the particle ends up with an orientation that is 
slightly tilted relative to the initial state. 

We next perform simulations on a channel twice as wide {W = 8D) and take four different values of 
Reynolds number. Re = 1.5, 2.4, 4.2, and 6.4. The particle trajectories are shown in Figure [T^ where we 
observe that in contrast with the W = AD results in Figure [TUl there is no overshoot of the centerline even 
for the highest value of Re. We attribute this behavior to the fact that in a wider channel, the hydrodynamic 
interactions between wall and particle that drive the horizontal motions are substantially weaker. The plot 
of angular velocity in Figure [T3] exhibits slightly different dynamics than the narrower channel, although the 
amplitude of the rotational motion is at least an order of magnitude smaller. This is to be expected since 
the rotational motion is also driven by the wall-particle interactions which are weaker for the wider channel 
case. 

We conclude our examination of the single-particle settling dynamics by comparing in Figure [14] the 
drag coefficients for the two different channel widths considered above, based on the formula Cd = t^{Ps — 
Pf)gD/{2pfV^) derived from equation (I5T|) . We have also included in this figure values of the drag coefficient 
taken from the following two papers: 



Feng, Hu and Joseph 16|, who performed numerical simulations using a finite element method for a 
single particle settling in channels of width AD and 8D. Because we will refer to this paper so often, 
we will refer to it with the abbreviation FHJ. 
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Figure 10: Settling trajectories for an off-center particle at different Reynolds number, in a channel of width W = AD. 



Sucker and Brauer [44J , who developed an empirical formula that is a fit to experimental data for the 
cylinder in a very large fiuid domain. They also developed an approximate analytical formula for an 
unbounded domain that matched closely with the experimental data. 



Our simulations match reasonably well with those of FHJ particularly for the W = 81? channel in the large 
Ly limit. Sucker and Brauer's empirical formula clearly deviates from both results because it applies strictly 
only to unbounded domains. While these results are encouraging, a much more comprehensive comparison 
is needed in order to draw any solid conclusions about the accuracy of our numerical approach. 



5. Numerical results: Two particle case. 

This section investigates the interactions between two circular particles with identical diameter D that 
are settling in a channel of length Ly = 3. We consider several initial configurations pictured in Figure [TSl 
and again compare the solution at different values of Reynolds number by varying viscosity over the range 
AtG [0.0015,0.16]. 

We begin by describing a well-known phenomenon in particle suspension flows wherein pairs of particles 
interact and undergo a "drafting, kissing and tumbling" behavior (which we abbreviate by DKT). This 
phenomenon has been established experimentally in papers such as 17|, |29| and demonstrated numerically 
in [lg|, and can be justified physically as follows. The leading particle creates in its wake a reduction of 
pressure as it falls under the influence of gravity. Provided that the trailing particle is close enough to 
interact with this wake, it experiences a smaller drag force than the leading particle. As a result the trailing 
particle falls faster and the particles approach each other - this is the initiation of the "drafting phase" . As 
the distance between the particles decreases, they eventually become close enough to nearly touch, which is 
referred to as "kissing" . The kissing particles momentarily form a single longer body that is aligned parallel 
with the flow; however, this parallel arrangement is unstable and the particles eventually tumble relative to 
each other and swap leading/trailing positions ~ this is the "tumbling phase" . The particles subsequently 
separate and one of two things happens: cither the DKT process repeats, or the particles continue to separate 
until the interaction force becomes so weak that they fall independently at their "natural" wall-corrected 
vertical settling velocity [38]. 

The simulations in this section are performed using the four initial configurations depicted in Figure [T51 

(a) Aligned vertically, one above the other along the channel centerline. 
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Figure 11: Plots of the angular velocity ixj in rad/s, for a single particle released from an off-center location, with W = 4D 
and Re = 4.9. Positive oj corresponds to counter-clockwise rotation, (a) Variation of angular velocity with time, (b) Angular 
velocity versus horizontal location, where the dashed line represents the channel centerline. 
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Figure 12: Settling trajectories for an off-center particle at different Reynolds number and in a channel of width W = 8D. 




Figure 13: Plot of angular velocity lu for a single particle released from an off-center location, with W = 8D and Re = 6.4. 
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Figure 14: Comparison of drag coefficients for a single particle settling in channels of width 4D and 8D, on a log— log scale. The 
experimental correlation of Sucker and Brauer |44| and numerical results from FHJ [lal are included for comparison purposes 



(reproduced with permission). Parameter values for our simulations: 
L^ = 0.68 (when W = 8D). 



D = 0.08, Ly = 3, Li, = 0.36 (when W = 4D) and 



(b) Aligned vertically, but shifted to the left to a position midway between the channel centerline and the 
left wall. 

(c) Aligned horizontally, and placed symmetrically about the centerline. 

(d) Aligned horizontally, but shifted to the left of center. 

In all cases, the particle diameter is D = 0.08 cm and the initial separation distance between the centers of 
mass of the two particles is 2D. As in the previous section, we also consider two channel widths, W — AD 
and W = 8D. 

5.1. Two vertically- aligned particles, released along the centerline. 

As a first test of the two-particle case, we use the initial set-up shown in Figure I15f a) wherein the 
particles are both released along the centerline with their centers of mass separated vertically by a distance 
2D. We perform simulations with channel width W = 81? and three different Reynolds numbers. Re = 3, 
14 and 80. 

Starting with the smallest value oi Re = 3, we find that both particles remain along the channel centerline 
throughout the simulation, and while drafting and kissing behavior is observed, no tumbling occurs. As seen 
in Figure [infa) , the trailing particle approaches quite close to the leading particle, but never touches it. 
This is because of the increase in the effective diameter of the particle owing to the delta function smoothing 
radius, as was discussed already at the end of section 14.3.21 After kissing, the particles continue to fall as 
a single body with no significant relative motion, except for a very slight "wobble" that corresponds to a 
small- amplitude oscillation in the orientation angle (refer to the angular velocity plot in Figure ITST b')'). This 
motion can be attributed primarily to oscillations of the IB points making up the particles that arises from 
the IB spring forces driving slight deviations from the equilibrium (stress-free) state; this is an artifact of 
the IB method that is not present in actual solid particles. 

Upon increasing the Reynolds number to Re — 14 the dynamics become much more complex. The 
trailing particle catches up with and subsequently passes the leading particle at time i « 11 s, breaking 
the left-right symmetry. This tumbling behavior is evident in Figure llTf a) where the vertical separation 
distance becomes negative. After the first tumble, the particles separate horizontally and move to locations 
symmetrically opposite each other relative to the centerline and separated by a horizontal distance of roughly 
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Figure 15: Four initial configurations for tfie two-particle settling problem, where the particle centers are separated by a 
distance 2D: (a) aligned vertically, along the channel centerlinc; (b) aligned vertically, offset to the left of center; (c) aligned 
horizontally, symmetric about the centerline; (d) aligned horizontally, offset to the left of center. 
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Figure 16: Settling dynamics of two particles initially aligned vertically, with parameters W = 8D and Re = 3. (a) Vertical 
separation distance, (b) Angular velocity ui (positive = counter-clockwise rotation). 
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0.27 cm. They continue to fall vertically at approximately the same x locations, and for the next 20 seconds 
they exchange leading and trailing positions via a small periodic variation in the vertical velocity whose 
amplitude decreases in time. By time i w 35 s, the particles have essentially reached a steady state in which 
they are falling at constant velocity and maintaining a constant separation distance (refer to Figure [ITjb)). 

The angular velocity plot in Figure [TTT c) shows that both particles experience a significant rotation 
during the tumbling phase that is several orders of magnitude larger than the small "wobbling" motion 
observed in the Re — 3 case. In fact, the growth of this rotational motion appears to be connected with the 
breaking of the horizontal symmetry that initiates the tumbling motion. By time i = 15 s, the rotational 
motion has subsided. 

Because of the symmetry in both the initial conditions and the governing equations, one would expect 
that the numerical solution should remain symmetric for all time, regardless of Reynolds number. The most 
likely source of asymmetry that initiates the tumbling behavior observed in the higher Re simulation is 
numerical error - these errors are sufficiently damped out when Re — 3, but remain large enough to initiate 
tumbling at Re = 14. This conjecture is borne out by the simulations in section [5^ where two particles are 
aligned vertically and released off-center. We have nonetheless shown the results for this symmetric initial 
condition since it is commonly studied in other simulations [16]. 

As the Reynolds number is increased yet further to Re = 80, we observe in Figure [TS] another qualitative 
change in solution behavior that is most easily seen in the sequence of snapshots collected in Figure [T21 The 
two particles begin with a DKT exchange such as that observed for Re = 14, however this occurs as the two 
particles drift together toward the left channel wall (instead of toward the channel center line). Following 
that, the particles drift reverse direction toward the right wall and undergo a second DKT exchange, after 
which the particle that was initially trailing ends up in the lead. These two DKT exchanges are accompanied 
by a back-and-forth rotational motion of each particle that has an amplitude similar in size to the Re = 14 
case (refer to Figure fTST b)). Once again, it appears to be growth in the small "wobble" in the particle 
angular velocity plot that initiates tumbling. We also observe that when a particle nears the left wall it 
experiences a clockwise rotation, while the direction of rotation is reversed near the right wall - this is 
consistent with physical intuition, which suggests that wall drag arising from the no-slip condition at the 
channel wall should cause a rolling-type motion as the portion of the particle closest to the wall slows down. 

FHJ \l&\ have performed a similar computation at Re = 70 (for the same symmetric initial conditions) 
that exhibits results consistent with ours up to time t « 5.4 s. However, they terminate their computation at 
this point and there is no indication in their paper of the subsequent dynamics. We have computed beyond 
this time and find that after the second tumble, the particles separate vertically as they migrate toward the 
center of the channel. By time t = 10 s they settle at roughly the same speed and no longer interact in 
any significant way . Our computations also exhibit much the same qualitative behavior as the experiments 
reported in IT], |29| and numerical simulations from 19|, |46[ . 



5.2. Two vertically- aligned particles, released off-center. 

In this section, we consider an asymmetric initial layout where the two particles are aligned vertically 
(again separated by a distance 2D) but instead have their centers of mass displaced to the left of center at 
location x — W/A. This initial geometry is depicted in Figure ITST b) . 

Results are first reported for a channel of width W = AD and four values of Reynolds number: Re — 1.5, 2, 
10 and 47. For the three smallest values of Reynolds number, we observe the behavior pictured in Figures [H] 
and l22l Initially, both particles drift toward the right, with the trailing particle moving toward the centerline 
while the leading particle moves to a point roughly mid-way between the centerline and the right wall. At 
the same time, the trailing particle speeds up in the wake of the leading particle so that they approach the 
same height. After this initial realignment, the two lowest Reynolds numbers [Re = 1.5, 2) undergo another 
slight adjustment in the horizontal locations so that the two particles are located symmetrically about the 
centerline; interestingly, the two particles in the i?e = 10 case remain in a slightly asymmetric layout relative 
to the centerline. After that, the particles fall with roughly constant speed and without changing a;-locations 
in the channel. The snapshots in Figure [22] for the Re = 2 case show that the particles do not enter either 
kissing or tumbling phases. 
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(b) Horizontal separation distance 
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Figure 17: Settling dynamics of two particles at Re = 14 initially aligned vertically along the centerline in a channel of 
width W = 8D. (a) Vertical separation distance, (b) Horizontal separation distance, (c) Angular velocity lj (positive = 
counter-clockwise) . 
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Figure 18: Settling dynamics of two particles at Re = 80 initially aligned vertically along the centerline in a channel of width 
W = 813. (a) Vertical separation distance, (b) Angular velocity ui (positive = counter-clockwise). 
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Figure 19: Snaphots of particle interactions in the channel of width W = 8D with Re ■ 
initially aligned vertically along the channel centerline. 
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Figure 20: Velocity vector plots for the drafting/kissing (left) and tumbling phase (right) for the case when the particles are 
initially aligned vertically along the channel centerline. Parameters: W = 8D and Re = 80. 

Both particles experience a distinct rotational motion as shown in Figure [53] for Re = 2, but the magni- 
tude of the angular velocity is not as large as was observed during the tumbling phase for the simulations 
in the previous section. Corresponding results for a channel of width W = 8D do not show any significant 
difference in qualitative behavior. 




Figure 21: Horizontal particle locations for two particles initially aligned vertically but off-center. Parameter values: W 
and Re = 1.5, 2, 10. 



4D 



A very different behavior is observed for the highest value of Reynolds number (Re = 47) as seen in 
Figures [M] and [55l Up to time i « 7 s the dynamics are similar to the lower Re cases in that the particles 
migrate to the right with the leading particle approaching closest to the wall. However, at this time the 
particles separate horizontally and a marked back-and-forth oscillation appears that grows in magnitude 
between times t « 7 and 23 s, until the oscillating particles overlap with each other near the centerline. 
At t « 23 s, the particles undergo a strong interaction in which they swap horizontal locations and move 
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Figure 22: Snapshots of particle interactions for the case when the particles are initially aligned vertically, but off-center. 
Parameter values: W = 4Z) and Re = 2. 
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Figure 23: Angular velocity w for two particles initially aligned vertically and off-center. Parameter values: W = AD and 
fie = 2. 
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to positions symmetrically opposite to each other on either side of the channel. At this stage, they have 
reached a steady state with roughly constant vertical velocity. The growing horizontal oscillations in the 
time interval [7, 23] are accompanied by a synchronized rotation of both particles (in opposite directions) 
that grows rapidly and then also dies out after t « 23 s (refer to Figure IMf b)). 

5.3. Two horizontally- aligned particles. 

We next simulate the motion of two particles initially aligned horizontally in a channel of width W — 8D. 
We consider two sets of initial conditions pictured in Figures [TSTc) and (d) , first where the particles are 
located symmetrically with respect to the channel centerline, and the second an asymmetric arrangement 
that is shifted to the left. 

Our main aim here is to determine to what extend our results are able to reproduce the finite element 
simulations of FHJ 16^ using W = 8D and Re = 1.52. They computed particle dynamics such as that 
pictured in Figure [211 that can be separated into three distinct phases: 

i. a first phase that consists of a rapid re-adjustment up to time t* « 500 (measured in dimensionless 

time units, with to t* = ty^g/D ) during which the particles separate horizontally to locations that 

are equally-spaced from the left and right walls, 
ii. a second phase where the particles maintain their horizontal positions and fall together with the same 

vertical speed until t* « 4000. 
iii. a third phase in which the particles shift together to the right into a new equilibrium state where the 

left-most particle oscillates about the centerline, while the right-most particle is much closer to the 

right wall and also oscillates side-to-side but with smaller amplitude. 

We remark that FHJ's simulations were intended to reproduce the experiments of Jayaweera and Mason [28[, 
wherein two long thin cylinders were settling in a large tank, with the same initial conditions and Re between 
0.1 and 1.0. Jayaweera and Mason's discussion of their experimental results makes mention of the first two 
phases but not phase iii. 

We begin with the symmetric case where the two particles have initial horizontal positions 

L^ W L W 

x = ^ 8" T T' ^^ 

and Reynolds number Re — 1.6 which is very close to FHJ's value. The horizontal locations of the simulated 
particles are pictured in Figure[27l with the plot axes rescaled to use the same dimensionless variables as FHJ 
in Figure [26l Our solution exhibits a steady state at long times that corresponds to positions x/D w 2.25 and 
5.75 located symmetrically across the centerline; these positions are very close to those obtained in FHJ's 
phase ii. Furthermore, the rapid transient in phase i ends at a dimensionless time of roughly t* — 500, 
which is also very close to FHJ's value. These two results suggest that our numerics are consistent with 
FHJ during phases i and ii and that we are capturing this portion of the motion properly. 

However, we do not capture the same phase iii behavior since our two particles never deviate from their 
steady state locations for t* > 4000 when FHJ's phase iii begins. Similar dynamics to FHJ's phase iii have 
also been computed by Aidun and Ding ^J with a lattice-Boltzmann method, and they ascribe this periodic 
behavior to the appearance of a solution bifurcation. It is likely that this bifurcation is sensitive not only to 
solution parameters but also to the presence of numerical error. Therefore, we suspect that the first order 
accuracy of our IB method may be preventing the numerics from capturing this transition to a periodic 
state at longer times. 

We repeated the previous calculation by increasing the Reynolds number to Re = 4.4 and our results are 
pictured in Figures B51 and [^ which exhibit similar dynamics to the lower Re case. In particular, we still 
observe no transition to phase iii behavior even at this higher Reynolds number. These results give us some 
confidence that our IB simulations are reproducing physically-relevant behavior corresponding to phases i 
and ii, but a more detailed numerical study is required in order to determine the source of the discrepancy 
between our method and FHJ's approach. 
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(a) Horizontal positions 




(b) Angular velocities 




(a) : leading particle 
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Figure 24; Settling dynamics of two particles initially aligned vertically and ofF-centcr, with parameters W = 4D and Re = 47. 
(a) Horizontal positions, (b) Angular velocity ui (positive = counter-clockwise). 
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Figure 25: Snapshots of particle interactions for two particles initially aligned vertically, but off-center. Read from top to 
bottom, then left to right. Parameter values: W = 4D and Re = 47. 
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Figure 26: Horizontal locations of the two particles, as simulated by FHJ. Parameters: W = SD, Re = 1.52. Reproduced from 
[la . Figure 32(b)], with permission. 
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Figure 27: Horizontal particle positions for two particles initially located on the same horizontal line, and symmetric about the 
centerline. Parameter values: W = SD and Re = 1.6. 
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Figure 28: Horizontal particle positions for two particles initially located on the same horizontal line. Both symmetrical and 
asymmetric initial conditions are pictured. Parameter values: W = 8D and Re = 4.4. 
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Figure 29: Snapshots of particle interactions for two particles initially aligned horizontally and centered. Parameter values: 
W = 8D and Re = 4.A. 
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The perfectly symmetric initial conditions used above are somewhat artificial, and will never actually 
occur in a real flow. Hence, we have also simulated an asymmetric initial placement of the particles given 

by 



X — — — and X = ——, 

2 4 2 ' 



(39) 



in which the initial particle locations from ([55| are shifted a distance W/^ to the left as pictured in Fig- 
urefTSld). Otherwise, the channel width W = 8D and Reynolds number Re = 4.4 remain the same as in the 
symmetric case. The numerical results are shown in Figures E51 and [501 where the particles exhibit similar 
dynamics to the symmetric case and approach the same long-term equilibrium solution. The only difference 
can be seen in the transient motion where the particles undergo one additional oscillation in the horizontal 
locations en route to steady state. 
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Figure 30: Snapshots of interactions for two particles initially aligned horizontally and off-center. Parameter values: W = 8D 
and Re = 4.4. 



6. Conclusions. 

The main aim of this paper is to demonstrate the ability of the immersed boundary method to simulate 
realistic dynamics of solid particles settling under gravity within a Newtonian incompressible fluid. The solid 
particles are modelled as a network of stiff springs, while the added mass of the particles is incorporated 
using an extra gravitational forcing term that is spread onto fluid points via a regularized delta function. 
Numerical simulations of a single particle show good agreement with the most accurate empirical formula for 
wall-corrected settling velocity due to Faxen. Furthermore, two-particle simulations reproduce qualitatively 
features of the dynamics seen in both experiments and numerical simulations. 

This study is by no means a comprehensive comparison to other results from the extensive literature 
on particle sedimentation, but rather sets the stage for such a study in future. In particular, we plan to 
perform a more detailed comparison with other published results, focusing first on our idealized cylindrical 
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particles. By implementing improvements to the numerical algorithm that increase accuracy of the solution 
approximation (such as in |.3ll|) we hope to be able to explain the discrepancy we observed between our 
results and those of Feng, Hu and Joseph [l6|. After that, the natural next step would be to extend our 
numerical method to 3D in order to permit simulations spherical particle interactions in a more realistic 
geometry. 

We emphasize that this study is a "proof-of-concept" that the immersed boundary method may be 
applied to simulating the sedimentation of particles that are denser than the suspending fluid. We make no 
claim to improve on or to compete with other numerical methods that are specially-tailored to deal with 
rigid, non-deformable particles. Instead, our ultimate goal is to solve sedimentation problems involving 
irregularly-shaped and highly deformable particles, which to our knowledge has not been sufficiently well 
studied in the literature. Such particle systems arise in the study of suspensions of red blood cells, wood pulp 
fibers, vesicles, bubbles, etc. Making use of the uniform triangulated meshes from the DistMesh package 
will allow us to deal with more general particle shapes. Furthermore, we plan to take advantage of recent 
developments in massively parallel immersed boundary algorithms by Wiens and Stockie |5l| , which should 
prove instrumental in allowing efficient 2D and 3D immersed boundary simulations to be performed for 
non-dilute suspensions containing large numbers of particles. 
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